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Fast symmetric matrix inversion using modified Gaussian elimination 


Anton Kochnev, Nikolai Savelov 


Abstract —In this paper we present two different variants 
of method for symmetric matrix inversion, based on modified 
Gaussian elimination. Both methods avoid computation of 
square roots and have a reduced machine time’s spending. 
Further, both of them can be used efficiently not only for 
positive (semi-) definite, but for any non-singular symmetric 
matrix inversion. We use simulation to verify results, which 
represented in this paper. 


I. INTRODUCTION 


Symmetric matrix inversion is one of the most important 
problem for many practical tasks e.g. analysis of electrical 
circuits with inductance elements [1], synthesis of Kalman 
or Wiener filters [2], using of finite element method [3]. 

Existing symmetric matrix inversion methods are 
Cholesky decomposition, LDL decomposition [4], bordering 
method [5], and the most efficient Krishnamoorthy - Menon’s 
method (based on Cholesky decomposition, and requires 
n 3 n 2 

-I-operations with n square roots computation) [6], 

[ 7 ]. 

The aim of this paper is to propose a symmetric matrix 
inversion method, which reduces machine time spending 
compared to Krishnamoorthy-Menon’s method by avoiding 
of square root computations. Moreover, this fact allow us to 
use the proposed method for efficient inverse of symmetric 
matrix not only with strict diagonal dominance, but without 
diagonal dominance as well. 


II. MODIFIED GAUSSIAN ELIMINATION 

In this section shows modified Gaussian elimination, 
which proposed method based on [8]-[12]. 

Let there is a system of linear equations (03, 


Let cti be an i th column of A; let F m be a matrix F after 
m th change; let f™ be an i th row of F m , m = 0, n; let 
F° = /, where I is an identity matrix. Then solution of © 
could be done with formulae © below: 

• /r +1 = //" f° r * < m +1, m = 1, n — 1, and x\ is an 

unrequired variable; 

. C = - fZ +1 for m = M = If 

fm+i a m+i = 0, then two rows of F should be 
permitted. Such permutation is always possible for non¬ 
singular matrix A (see [8]-[ 12]); 

• /r +1 = fr - (fr a rn+i)f™ti for other re q uired *; 

• Xi = ffb for any required Xi. 

( 2 ) 


It can easily be checked that /"a, = 1, /"a,- = 0, for 
any required Xi and for j = 1 ,n, j =4 i. Further, if all 
elements of x are required variables, then F" = A ~ 1 . Ma¬ 
trix inversion using modified Gaussian elimination requires 
n 3 multiplications and divisions. If only the last element 
x n of x is a required variable, then Gaussian elimination 

Tl 3 Tl Tl 

requires-1-1— multiplies and divisions. Let p be a 

3 2 6 _ 

number of required variables in x: Xi, i = n — p + 1, n is a 

required variable. Then number of multiplies and divisions 

for Gaussian elimination could be determined with formula 

n 3 nr n 2 p 3 p 2 p , 

T + T + 6 +pn - pn ~T + Y-6 (see [12]) ' 

III. PROPOSED METHOD 

In this section shows a method for efficient symmetric 
matrix inversion based on modified Gaussian elimination. 


Ax = b, (1) 

where A £ C nxn , x £ C raxl , b £ C nxl . During modified 
Gaussian elimination an addition matrix F : F £ C™ x " 
changes instead of matrix A, but addition memory is not 
necessary in this case [12], 

Let vector x consist of two types of variables: required, 
which should be find during elimination, and unrequired, 
which is not interesting for researcher. Re-solution of © 
after some changes of A could be done with reduced number 
of multiplications and divisions, using formulae from [8]- 
[ 12 ]. 
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A. First variant of the proposed method 

The first variant of the proposed method consist of two 
different stages. On the first stage we use formulae (2) for 
p = 1, and only x n is a required variable (indeed, it is a 
valid proposition for p = 0 as well). On the second stage we 
use addition formulae described below. 

Let us introduce some notation. 

Let /™ +1 be an element from i th row and j ih column 
of F m+l \ let ctij be an element from i th row and j th 
column of A; let Sy +1 be a submatrix of F m+1 , such that 
sy +1 = (/" l+1 ), * = l,?7i + l, j = l,m+l; let S a be 
a submatrix of A, such that S a = {ctij), i = l,m+l, 

j = 1, 777 + 1. 

It is easily shown that after the first stage F is a lower 
triangular matrix. After the second stage F became a matrix 
A -1 , as if we use (2) for p = n. 










For symmetric matrix inversion using modified Gaussian 
elimination that is enough to use only lower triangular 
matrix. To prove this statement, we need a lemma [Q 

Lemma 1: Let F m+1 is a matrix from © for p = n, 

m = 0,n- 1. Then SJ+ 1 = (S'a) -1 - 
Proof: 

Using 0 , let us consider 4 cases. 

Case 1. Let i = m + 1. Then 

fm+1 _ rm+1 

Ji u i — J m+ l u m+l) 

fm+1 _ ___ fm _ 1 

Jm+l u m+l — r m Jm+l u m+l — r, 

Jm+l u m+l 

m = 0, n — 1; 

Case 2. Let j = m + 1. Then 

fm+1 _ fm+1 
Ji a j ~ Ji a m+ 1, 

fm+1 _ f m „ fm fm+1 _ O 

Ji a m+l Ji &m +1 Ji (lm+lj m -fi &m+l — U, 

m = 0, n — 1, i = 1, m; 


Case 3. Let i = to. Then 
f m+1 a — f m+1 n 

Ji u i Jm u m, 

fm+1 _ fm _ fm fm+1 

Jm u m~Jm u m 


fmUm+1 


fm+ L „ fm „ jrn^m-r l fm „ i. 

Jm a m ~ Jm a m J m+l a m ~ 

Jm+ l a m+l 


It can be checked easily for Vi : * < to + 1, m = 0, n — 1; 

Case 4. Let j = m. Then 

f m +l _ f m n _ frn f m +l n _ n. 

Ji u j Ji u m Ji Um+lJ m+l a m — O, 

It can be checked easily for Vj: j < m + 1, j ^ i, 
to = 0, n — 1. 

■ 

Since lemma [j] it follows that S'™ +1 is a symmetric 
matrix. 

This statement is needed for the second stage of the 
method. Let F k be a matrix F after k th change, k = 0, n, 
where F° is a lower triangular matrix after the first stage, 
and F n = A~ x . Let f\ k is an element from i th row and 
j th column of F k . Then the second stage describes with 
formulae below: 


fij = fij ^or i = k + 1, n, j = 1 , n, k = 0, n - 1, 
and for i = 1, k, j > i, k = 1, n — 1; 

rk+l rk+1 

f-J 1 = '— for i = lfk, j = M, 


f, 


fc+i 

fc+ifc+i 


k = 1, n — 1. 


( 3 ) 

It is easily to prove that A _1 = F n + ( F n — D) T , where 
D is a diagonal matrix, such that da = i = l,n. 

Number of multiplications and divisions for the first stage 

1 9 

n n n 

of the method describes with formula-1-1—, as we 

3 2 6 

note earlier. 

It is easily shown that number of multiplications and 
divisions for the second stage describes with formula 

n 3 n 2 2 n 

IT + ~2 3~ 


Both stages for symmetric matrix inversion using the 

n 3 n TL 

first variant of the proposed method requires — + n — — 
multiplications and divisions. It is less then requirements of 
Cholesky decomposition or LDL decomposition (see table 
1 ). 

Let us remark that proposed method avoid square root 
computations; this considerably reduce machine time spend¬ 
ing, and make it possible to use proposed method not only 
for positive determined, but for any invertable symmetric 
matrices as well. 

B. Second variant of the proposed method 

The second variant of the proposed method consist of only 
one stage. 

Suppose, that 


pm+l = pm + AF m ; (4) 

where F™ = F m , but all elements from (to + l) ih row 
of FJ" are zeros. 

Let /"' be an element of i th row and j th column of F rn ; 
let fl'ij be a j th column of F m . 

For explanation of the following formulae, we need a 
lemma 

Lemma 2: Let F m+1 is a matrix from (0 for p = n, m = 
0, n — 1. Then A F m from 0 be A F m = 

Proof: From © A F m = F m+1 - F'f. 

Using 0, we get 

1. 

fm+l _ _^_ fm 

J m+1 rm J m-\- 1’ 

J m+lUrn+l 

rm -\-1 fm -\-1 rm 

J m+l J m+lm+lJ m+li 

fmtlm+lfm+l = lmt\ ~ 0, m = 0^1; 

2 . 


fm+1 


rm+1 rm / rm \ rr 

Ji Ji \J i 

rm+1 _ rm 1 rm 

Ji Ji fm n /m+l 5 

Jm+l a m+l 

rm+1 rm / fm \ rm rm 

Ji Ji \Ji tt m+lJ/m-|-lm+l/m+l’ 

fm+1 _ rm . fm+1 rm 

Ji Ji ' Jim+lJm+li 

_ f m 

J i ’ 


fm+1 rm _ fm+1 

J im+lJ m+l Ji 


to = 0, n — l,i = l,n, i m + l. 


If we combine lemma Q] with lemma 01 we get formulae 
below: 


f 

J 1 


m+l 

im+1 


1 


fm n 

J m+l^m+l 


/; 


m+li 


for to = 0, n — 1, 


i = 1 , to + 1 , 

fZ + +\ = - (fram+i)fZtL+i for rn = 0 ^ 2 , 

i = to + 2, n, 

f™ +1 = fij+flUJifm+ij for m = l,n- 1, j = l~m, 
i = j,n, i to + 1, 

fZXlj = fjm+i for m = I^T, j = l~m. 





























• /™ +1 = f™ +1 for m = 2, n — 1, j = 2 ,m, 

* = _ _ _ 

• fij +l = flj for m = 0, n — 2, i = 1, n, j = m + 2, n. 

(5) 

The formulae (0 describe an idea of the method in detail, 
but for practical tasks it is better to use different formulae: 

• fZtlj = 7m ~- fm+n for m = 0,n-l, 

_ Jm +1 a m +1 

j = 1, m + 1, 

• /im +1 = ~ (/r tt m+l)/m+lm+l f° r m = 0,71—2, 
i = to + 2, n, 

• /^ +1 = /y + fZWfZ+i! for m = 1^2, i = 
to + 2, n, j = 1, to, 

• fiJ+Ctufm+ij for m = T^T, i = 1/m, 
j = M> 

• /" +1 = f ™ +1 for to = 0,n-2, * = l,n, 
j = to + 2, n, 

• fij +1 = fij for TO = 1,71— 1, * = 1, TO, 

j = 7 + 1, TO + 1. 

(6) 

It can easily be shown that 0 and 0 are equivalent. 

Number of multiplications and divisions for the second 
variant of the proposed method describes with formulae 

n 3 n 2 

— + —. It is less then requirements of Cholesky decom¬ 
position, LDL decomposition or Krishnamoorthy - Mennon 
method. It is the same requirements as for bordering method, 
but it should be noted that bordering method could not be 
use for inversion of matrix with Mu ^ 0, i = l,n, where 
Mij is a minor of A, i = 1, n, j = 1, n. 

At the same time, proposed method could be use for 
inversion of matrix with Mu = 0, if A is not a singular 
matrix. 

Let us remark that the second variant of the proposed 
method avoid square root computations as well. 

IV. SIMULATION SETUP 

In order to demonstrate advantages of the proposed al¬ 
gorithms, we use MATLAB based simulation via different 
CPUs. We give results for Intel Core i5-3230M 2.60 GHz 
below (Intel Pentium Dual Core T 2390 1.86 GHz and 
Intel Atom N450 1.67 GHz gives familiar results). We 
compare proposed algorithms with the most efficient notable 
algorithms for symmetric matrix inversion: Cholesky decom¬ 
position [13], LDL decomposition [4], and Krishnamoorthy- 
Mennon method [6], [7], [14], We generate table with full 
equations, which describes number of multiplications, divi¬ 
sions and square roots computation for every noted method. 

The first row of the table 1 describes Cholesky decomposi¬ 
tion and solving of systems of linear equations LB = J, and 
L t A~ 1 = B. The second row describes LDL decomposition 
and solving of SLE LX = /, DB = X, and L T A~ 1 = B. 
The third row describes matrix inversion using Krishnamoor¬ 
thy -Mennon method, based on Cholesky decomposition [6], 
[7]. The fourth row describes the first variant, and the fifth 
row describes the second variant of the proposed method. 


Experiment 1. Inversion of a real symmetric matrix with 
strict diagonal dominance. Let qtheor be a number of mul¬ 
tiplications and divisions, and s t heor be a number of square 
root computations, determined with formulae from the table 
1. Let q pr act and s prac t be numbers of operations, determined 
with counter variables from MATLAB scripts. Results of the 
experiment are given in tables 2 and 3. 

Experiment 2. Inversion of real symmetric matrices of 
order n with strict diagonal dominance. Let f be a time for 
matrix inversion; let norm be a second norm : norm = 

11 A” 1 — Au^ v 112, where A” 1 is a matrix, inverted via one 
of described methods, and A~ 3 V is a matrix, inverted via 
MATLAB function inv (A). Results of the experiment are 
given in tables 4 and 5. 

Experiment 3. Inversion of real symmetric matrices of or¬ 
der n without diagonal dominance. Results of the experiment 
are given in the table 6. 

V. SIMULATION RESULTS 

From tables 2 and 3 we can conclude that formulae from 
the table 1 are correct. 

From tables 4 and 5 we can conclude that both variants of 
the proposed method provide notable reduction of machine 
time spending and has a good accuracy. 

From the table 6 we can conclude that both variants of the 
proposed method increase advantages for matrices without 
diagonal dominance. Let us remark that it is especially 
important for inductance matrix inversion. 

VI. CONCLUSIONS 

We propose a new method for symmetric matrix inversion 
based on modified Gaussian elimination with avoiding of 
square root computations. Proposed method could be useful 
for any scientific and technical problem with symmetric 
matrix inversion, especially if matrix has not a diagonal 
dominance. 






















METHOD OF MA¬ 
TRIX INVERSION 

NUMBER OF 

MULTIPLICATIONS 
AND DIVISIONS 

NUMBER OF 

SQUARE ROOT 

COMPUTATIONS 

Cholesky decomposi¬ 
tion 


n 

LDL decomposition 

SOiiliiB 

0 

Krishnamoorthy 
Mennon’s method 

(based on Cholesky 
decomposition) 

ri 6 n 2 

T + T 

n 

The first variant of the 
proposed method 

IQS 

0 

The second variant of 
the proposed method 


0 


Table 1. Table summarizing the number of operations for inversion of a 
matrix with strict diagonal dominance via different methods. 


METHOD OF MA¬ 
TRIX INVERSION 

Qtheor 

Qpractr 

$theor 

Spract 

Cholesky decomposi¬ 
tion 

515000 

515000 

100 

100 

LDL decomposition 

671650 

671650 

0 

0 

Krishnamoorthy 
Mennon’s method 

(based on Cholesky 
decomposition) 

505000 

505000 

100 

100 

The first variant of the 
proposed method 

509950 

509950 

0 

0 

The second variant of 
the proposed method 

505000 

505000 

0 

0 


Table 2. Table summarizing the number of operations for inversion of a 
matrix with strict diagonal dominance via different methods (n = 100). 


METHOD OF MA¬ 
TRIX INVERSION 

Qtheor 

Qpractr 


Spract 


62875000 

62875000 

500 

500 

LDL decomposition 

83458250 

83458250 

0 

0 

Krishnamoorthy 
Mennon’s method 

(based on Cholesky 
decomposition) 

62625000 

62625000 

500 

500 

The first variant of the 
proposed method 

62749750 

62749750 

0 

0 

The second variant of 
the proposed method 

62625000 

62625000 

0 

0 


Table 3. Table summarizing the number of operations for inversion of a 
symmetric matrix with strict diagonal dominance via different methods (n = 
500). 


METHOD OF MA¬ 
TRIX INVERSION 

n = 

100 

n = 

300 

n = 

500 

n = 

1000 

Based on Cholesky 
decomposition 

0.384 s. 

9.499 s. 

41.96 s. 

329.7 s. 

LDL decomposition 

0.332 s. 

7.950 s. 

33.23 s. 

261.0 s. 

Krishnamoorthy 
Mennon’s method 

(program by A. 
Krishnamoorthy [14]) 

0.125 s. 

2.950 s. 

12.12 s. 

96.84 s. 

Krishnamoorthy 
Mennon’s method 

(program with 

element-by-element 
access) 

0.113 s. 

2.793 s. 

12.44 s. 

103.3 s. 

The first variant of the 
proposed method 

0.046 s. 

0.743 s. 

3.086 s. 

27.13 s. 

The second variant of 
the proposed method 

0.046 s. 

0.718 s. 

2.845 s. 

25.64 s. 


Table 4. Table summarizing times of numerical computations for inversion 
of a symmetric matrix with strict diagonal dominance via different methods. 


METHOD OF MA¬ 
TRIX INVERSION 

n = 

100 

n = 

300 

n = 

500 

n = 

1000 

Based on Cholesky 
decomposition 

9.4E-19 

1.3E-18 

1.9E-18 

3.0E-18 

LDL decomposition 

1.0E-18 

1.5E-18 

2.0E-18 

3.0E-18 

Krishnamoorthy 
Mennon’s method 

(program by A. 
Krishnamoorthy [14]) 

9.4E-19 

1.3E-18 

1.9E-18 

3.0E-18 

Krishnamoorthy 
Mennon’s method 

(program with 

element-by-element 
access) 

9.4E-19 

1.3E-18 

1.9E-18 

3.0E-18 

The first variant of the 
proposed method 

1.5E-18 

2.5E-18 

4.1E-18 

6.7E-18 

The second variant of 
the proposed method 

1.5E-18 

2.5E-18 

4.1E-18 

6.7E-18 


Table 5. Table summarizing || — A—^ v \\ 2 of numerical computations 

for inversion of a symmetric matrix with strict diagonal dominance via 
different methods. 


METHOD OF MA¬ 
TRIX INVERSION 

n = 

100 

n = 

300 

n = 

500 

n = 

1000 

Based on Cholesky 
decomposition 

0.448 s. 

11.16 s. 

50.47 s. 

401.2 s. 

LDL decomposition 

0.328 s. 

7.420 s. 

33.22 s. 

264.4 s. 

Krishnamoorthy 
Mennon’s method 

(program by A. 
Krishnamoorthy [14]) 

0.153 s. 

6.185 s. 

54.42 s. 

861.7 s. 

Krishnamoorthy 
Mennon’s method 

(program with 

element-by-element 
access) 

0.318 s. 

10.66 s. 

73.34 s. 

1021 s. 

The first variant of the 
proposed method 

0.045 s. 

0.742 s. 

3.093 s. 

27.05 s. 

The second variant of 
the proposed method 

0.045 s. 

0.701 s. 

2.836 s. 

25.30 s. 


Table 6. Table summarizing times for inversion of a symmetric matrix 
without diagonal dominance via different methods. 
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